Generalidades

La calibración del modelo se hace de forma conjunta en términos de cantidad y calidad de agua utilizando una función objetivo que indica el grado de similitud entre los valores simulados y observados de caudal diario (m³/s), la concentración diaria de nitrato (NO3) y la concentración diaria de fósforo total (PT). Para establecer las concentraciones diarias de PT y NO3 se utilizaron las relaciones exponenciales que se presentan en la sección Sistema de Monitoreo:muestreo de fósforo_y_nitrato

Para ello se utilizaron los siguientes paquetes de R: hydroPSO y SWATplusR.

El código utilizado para la calibración se muestra en el Anexo.

Función objetivo

La función objetivo utiliza la eficiencia de Kling-Gupta (KGE). Se aplica a: caudal, fósforo y nitrato. Luego se añade una eficiencia relativa al sesgo de fósforo y nitrato. Finalmente se aplican pesos para dar mas importancia a el caudal así como a la parte baja de la cuenca. Esto se hace por dos razones:

• Primero: para suviazar los errores de observación que suelen ocurrir en cuencas pequeñas, estos errores son producto de la mayor variabilidad del sistema hidrológico a dichas escalas.
• Segundo: para dar mayor relevancia a la observacion diaria (en el caso de caudal) y disminuir el efecto que pudiese tener la regresion de fósforo y nitrato que ha sido utilizada.

La funcion objetivo (FO) se define por:

F_obj = 0.4F_obj|_x = 23 + 0.6F_obj|_x = 41

Donde 0.4 y 0.6 son los pesos que se aplican a las funciones objetivo (F_obj) en cada punto de control (x). Para Paso de los Troncos x = 23 y Paso Pache x = 41.

F_obj|_x se define de la siguiente manera:

F_obj|_x = 0.5KGE(Flow)|_x + 0.1KGE(PT)|_x + 0.1KGE(NO3)|_x + 0.15PBE(PT)|_x + 0.15PBE(NO3)|_x

Donde KGE es la eficiencia de Kling-Gupta, PBE es una eficiencia relativa al sesgo y los coeficientes 0.5, 0.1 y 0.15 son los pesos utilizados.

$$KGE = 1- \sqrt{(r-1)^2 + (\frac{\sigma_{sim}}{\sigma_{obs}}-1)^2 + (\frac{\mu_{sim}}{\mu_{obs}}-1)^2}$$

r es la correlacion lineal entre las observaciones y las simulaciones. σ_obs es la desviación estándar en las observaciones. σ_sim la desviación estándar en las simulaciones. μ_obs es la media de las observaciones, μ_sim el promedio de las simulaciones. Esto fue calculado utilizando la función KGE del paquete hydroGOF.

Luego, PBE se define como uno menos el valor absoluto de la suma de los residuos dividida entre la suma de las observaciones.

$$PBE = 1 - abs\left(\frac{\sum(sim-obs)}{\sum(obs)}\right)$$

Parámetros a optimizar

Se consideraron 22 parámetros en cada region. La cuanca alta está definida por las subcuencas “1,8,10,14,15,16,17,18,19,20,22” y la cuanca baja por “2,3,4,5,6,7,9,11,12,13,21,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41”. El rango de variación de cada parámetro así como el tipo de cambio se muestra en la siguiente tabla:

par_all <- readRDS("~/R_projects/SWAT_calibration_ETmodis/resultados/optim20200401_20/par_all.RDS")
best_sim <- readRDS("~/R_projects/SWAT_calibration_ETmodis/resultados/optim20200401_20/best_sim.RDS")

par_upr = apply(par_all,2,max)
par_lwr = apply(par_all,2,min)

par_upr = par_upr[best_sim$parameter$definition$par_name]
par_lwr = par_lwr[best_sim$parameter$definition$par_name]

par_optim = data.frame(parametro=best_sim$parameter$definition$par_name,
                       archivo=best_sim$parameter$definition$file_name,
                       cambio=best_sim$parameter$definition$change,
                       L_inf = as.numeric(round(par_lwr,2)),
                       L_sup = as.numeric(round(par_upr,2)),
                       Best_par = as.numeric(round(best_sim$parameter$values,2)))


kable(par_optim, caption ="Tipo de cambio, rangos de los parámetros en la optimización y mejor parámetro") %>%
  kable_styling(c("striped", "bordered"))

Resultados de la optimización

best_skill <- readRDS("~/R_projects/SWAT_calibration_ETmodis/resultados/optim20200401_20/best_skill.RDS")
library(plotly)
fig <- plot_ly(y = ~best_skill[,"KGE41"], type = "box",name="Paso Pache \n(calibración)" )
fig <- fig %>% add_trace(y = ~best_skill[,"KGE36"],name="Fray Marcos \n(validación)")
fig <- fig %>% add_trace(y = ~best_skill[,"KGE23"],name="Paso de los Troncos \n(calibración)")
fig <- fig %>% layout(yaxis = list(title="KGE [-]"), showlegend = FALSE)
fig

fig <- plot_ly(y = ~best_skill[,"pbias41"], type = "box",name="Paso Pache \n(calibración)" )
fig <- fig %>% add_trace(y = ~best_skill[,"pbias36"],name="Fray Marcos \n(validación)")
fig <- fig %>% add_trace(y = ~best_skill[,"pbias23"],name="Paso de los Troncos \n(calibración)")
fig <- fig %>% layout(yaxis = list(title="BIAS [%]"), showlegend = FALSE)
fig

library(zoo)
obs.data <- readRDS("~/R_projects/SWAT_calibration_ETmodis/data/obs.data.RDS")
obs.data = data.frame(obs.data)
obs.data$q23lag1 = c(obs.data$q23[2:nrow(obs.data)],NA)
obs.data$q36lag1 = c(obs.data$q36[2:nrow(obs.data)],NA)
obs.data$q41lag2 = c(obs.data$q41[3:nrow(obs.data)],NA,NA)

Qobs = obs.data$q23lag1
Qsim = best_sim$simulation$FLOW_23
Rdate = best_sim$simulation$date

fig <- plot_ly(x = ~Rdate, y = ~Qobs, type="scatter",mode = "lines",
               name="Qobs", line = list(color = 'black', width = 2))

fig <- fig %>% add_lines(x = ~Rdate, y = ~Qsim, name="Qsim",
                         type = 'scatter',
                         mode = 'line', 
                         line = list(color = 'orange', width = 2))

fig <- fig %>% layout(legend = list(x = 0.1, y = 0.9),yaxis = list(title="Flow (m3/s)"),xaxis = list(title="") )
fig

Qobs = obs.data$q36lag1
Qsim = best_sim$simulation$FLOW_36
Rdate = best_sim$simulation$date

fig <- plot_ly(x = ~Rdate, y = ~Qobs, type="scatter",mode = "lines",
               name="Qobs", line = list(color = 'black', width = 2))

fig <- fig %>% add_lines(x = ~Rdate, y = ~Qsim, name="Qsim",
                         type = 'scatter',
                         mode = 'line', 
                         line = list(color = 'orange', width = 2))

fig <- fig %>% layout(legend = list(x = 0.1, y = 0.9),yaxis = list(title="Flow (m3/s)"),xaxis = list(title="") )
fig

Qobs = obs.data$q41lag2
Qsim = best_sim$simulation$FLOW_41
Rdate = best_sim$simulation$date

fig <- plot_ly(x = ~Rdate, y = ~Qobs, type="scatter",mode = "lines",
               name="Qobs", line = list(color = 'black', width = 2))

fig <- fig %>% add_lines(x = ~Rdate, y = ~Qsim, name="Qsim",
                         type = 'scatter',
                         mode = 'line', 
                         line = list(color = 'orange', width = 2))

fig <- fig %>% layout(legend = list(x = 0.1, y = 0.9),yaxis = list(title="Flow (m3/s)"),xaxis = list(title="") )
fig

fig <- plot_ly(y = ~best_skill[,"kge_PT_41"], type = "box",name="Paso Pache \n(calibración)" )
fig <- fig %>% add_trace(y = ~best_skill[,"kge_PT_23"],name="Paso de los Troncos \n(calibración)")
fig <- fig %>% layout(yaxis = list(title="KGE [-]"), showlegend = FALSE)
fig

fig <- plot_ly(y = ~best_skill[,"pbias_PT_41"], type = "box",name="Paso Pache \n(calibración)" )
fig <- fig %>% add_trace(y = ~best_skill[,"pbias_PT_23"],name="Paso de los Troncos \n(calibración)")
fig <- fig %>% layout(yaxis = list(title="Bias [%]"), showlegend = FALSE)
fig

# datos observados
area_41 = 4896
area_36 = 2744
area_23 = 687
qlt_41 <- readRDS("~/R_projects/SWAT_calibration_ETmodis/data/rch41_calidad.RDS")
qlt_23 <- readRDS("~/R_projects/SWAT_calibration_ETmodis/data/rch23_calidad.RDS")
qlt_41$qobs = 1000*qlt_41$qobs/area_41
qlt_23$qobs = 1000*qlt_23$qobs/area_23
qlt_41$NO3_mgxs = qlt_41$NO3_mgxs/area_41
qlt_23$NO3_mgxs = qlt_23$NO3_mgxs/area_23
qlt_41$PT_mgxs = qlt_41$PT_mgxs/area_41
qlt_23$PT_mgxs = qlt_23$PT_mgxs/area_23
fit41_PT = lm(log(qlt_41$PT_mgxs)~log(qlt_41$qobs))
fit41_NO3 = lm(log(qlt_41$NO3_mgxs)~log(qlt_41$qobs))
fit23_PT = lm(log(qlt_23$PT_mgxs)~log(qlt_23$qobs))
fit23_NO3 = lm(log(qlt_23$NO3_mgxs)~log(qlt_23$qobs))
fit41_PT = round(c(exp(fit41_PT$coefficients[1]), fit41_PT$coefficients[2]),2)
fit41_NO3 = round(c(exp(fit41_NO3$coefficients[1]), fit41_NO3$coefficients[2]),2)
fit23_PT = round(c(exp(fit23_PT$coefficients[1]), fit23_PT$coefficients[2]),2)
fit23_NO3 = round(c(exp(fit23_NO3$coefficients[1]), fit23_NO3$coefficients[2]),2)
qlt_41$fit_PT = fit41_PT[1]*qlt_41$qobs^fit41_PT[2]
qlt_41$fit_NO3 = fit41_NO3[1]*qlt_41$qobs^fit41_NO3[2]
qlt_23$fit_PT = fit23_PT[1]*qlt_23$qobs^fit23_PT[2]
qlt_23$fit_NO3 = fit23_NO3[1]*qlt_23$qobs^fit23_NO3[2]

# datos simulados
qsim = best_sim$simulation

qlt.sim_41 = data.frame(flow_lps_km2=1000*qsim$FLOW_41/4896,
                        PT_mgs_km2 = 1000000*qsim$PT_41/(24*60*60*4896),
                        NO3_mgs_km2 = 1000000*qsim$NO3_41/(24*60*60*4896),
                        qobs_lps_km2 = 1000*obs.data$q41lag2/4896,
                        PTobs_mgs_km2 = 0.18*(1000*obs.data$q41lag2/4896)^1.02,
                        NO3obs_mgs_km2 = 0.11*(1000*obs.data$q41lag2/4896)^1.41)

qlt.sim_41[qlt.sim_41<=0] = NA
qlt.sim_41 = na.omit(qlt.sim_41)

qlt.sim_23 = data.frame(flow_lps_km2=1000*qsim$FLOW_23/687,
                    PT_mgs_km2 = 1000000*qsim$PT_23/(24*60*60*687),
                    NO3_mgs_km2 = 1000000*qsim$NO3_23/(24*60*60*687),
                    qobs_lps_km2 = 1000*obs.data$q23lag1/687,
                    PTobs_mgs_km2 = 0.1*(1000*obs.data$q23lag1/687)^0.93,
                    NO3obs_mgs_km2 = 0.08*(1000*obs.data$q23lag1/687)^1.41)

qlt.sim_23[qlt.sim_23<=0] = NA
qlt.sim_23 = na.omit(qlt.sim_23)

fit = lm(log(PT_mgs_km2)~log(flow_lps_km2), data=qlt.sim_41)
a = exp(fit$coefficients[1])
b = fit$coefficients[2]

q41fit = seq(min(qlt.sim_41$flow_lps_km2), max(qlt.sim_41$flow_lps_km2), length.out = 30)
PT41fit = a*q41fit^b
 

PT = plot_ly(data = qlt.sim_41, x=~flow_lps_km2,y=~PT_mgs_km2, type='scatter', 
             mode='markers', name="Simulado",
             marker=list(color='gray', size=4))

PT = PT %>% add_markers(data=qlt_41,x = ~qobs, y = ~PT_mgxs,
                        type = 'scatter',
                        mode = 'markers', 
                        marker = list(color = 'red', size=8),
                        name="Observado")

PT <- PT %>% add_lines(data=qlt_41,x = ~qobs, y = ~fit_PT,
                       type = 'scatter',
                       mode = 'line', 
                       marker = list(color = 'red', size=0),
                       line = list(color = 'red', width = 2),
                       name="Regresión observados")

PT <- PT %>% add_lines(x= ~q41fit, y = ~PT41fit,
                       type = 'scatter',
                       mode = 'line', 
                       marker = list(color = 'black', size=0),
                       line = list(color = 'black', width = 2),
                       name="Regresión simulados")

PT <- layout(PT, xaxis = list(title="Caudal (lps/km2)",range=c(-1,2),type = "log"),
             yaxis = list(title="PT (mg/s/km2)", range=c(-2,1),type = "log"),
             legend = list(x = 0.01, y = 1)
             #               showlegend = FALSE,
)


PT

fit = lm(log(PT_mgs_km2)~log(flow_lps_km2), data=qlt.sim_23)
a = exp(fit$coefficients[1])
b = fit$coefficients[2]

q23fit = seq(min(qlt.sim_23$flow_lps_km2), max(qlt.sim_23$flow_lps_km2), length.out = 30)
PT23fit = a*q23fit^b


PT = plot_ly(data = qlt.sim_23, x=~flow_lps_km2,y=~PT_mgs_km2, type='scatter', 
             mode='markers', name="Simulado",
             marker=list(color='gray', size=4))

PT = PT %>% add_markers(data=qlt_23,x = ~qobs, y = ~PT_mgxs,
                        type = 'scatter',
                        mode = 'markers', 
                        marker = list(color = 'red', size=8),
                        name="Observado")

PT <- PT %>% add_lines(data=qlt_23,x = ~qobs, y = ~fit_PT,
                       type = 'scatter',
                       mode = 'line', 
                       marker = list(color = 'red', size=0),
                       line = list(color = 'red', width = 2),
                       name="Regresión observados")

PT <- PT %>% add_lines(x= ~q23fit, y = ~PT23fit,
                       type = 'scatter',
                       mode = 'line', 
                       marker = list(color = 'black', size=1),
                       line = list(color = 'black', width = 2),
                       name="Regresión simulados")

PT <- layout(PT, xaxis = list(title="Caudal (lps/km2)",range=c(-1,2),type = "log"),
             yaxis = list(title="PT (mg/s/km2)", range=c(-2,1),type = "log"),
             legend = list(x = 0.01, y = 1)
             #               showlegend = FALSE,
)


PT

fig <- plot_ly(y = ~best_skill[,"kge_NO3_41"], type = "box",name="Paso Pache \n(calibración)" )
fig <- fig %>% add_trace(y = ~best_skill[,"kge_NO3_23"],name="Paso de los Troncos \n(calibración)")
fig <- fig %>% layout(yaxis = list(title="KGE [-]"), showlegend = FALSE)
fig

fig <- plot_ly(y = ~best_skill[,"pbias_NO3_41"], type = "box",name="Paso Pache \n(calibración)" )
fig <- fig %>% add_trace(y = ~best_skill[,"pbias_NO3_23"],name="Paso de los Troncos \n(calibración)")
fig <- fig %>% layout(yaxis = list(title="Bias [%]"), showlegend = FALSE)
fig

fit = lm(log(NO3_mgs_km2)~log(flow_lps_km2), data=qlt.sim_41)
a = exp(fit$coefficients[1])
b = fit$coefficients[2]

q41fit = seq(min(qlt.sim_41$flow_lps_km2), max(qlt.sim_41$flow_lps_km2), length.out = 30)
NO3_41fit = a*q41fit^b


PT = plot_ly(data = qlt.sim_41, x=~flow_lps_km2,y=~NO3_mgs_km2, type='scatter', 
             mode='markers', name="Simulado",
             marker=list(color='gray', size=4))

PT = PT %>% add_markers(data=qlt_41,x = ~qobs, y = ~NO3_mgxs,
                        type = 'scatter',
                        mode = 'markers', 
                        marker = list(color = 'red', size=8),
                        name="Observado")

PT <- PT %>% add_lines(data=qlt_41,x = ~qobs, y = ~fit_NO3,
                       type = 'scatter',
                       mode = 'line', 
                       marker = list(color = 'red', size=0),
                       line = list(color = 'red', width = 2),
                       name="Regresión observados")

PT <- PT %>% add_lines(x= ~q41fit, y = ~NO3_41fit,
                       type = 'scatter',
                       mode = 'line', 
                       marker = list(color = 'black', size=1),
                       line = list(color = 'black', width = 2),
                       name="Regresión simulados")

PT <- layout(PT, xaxis = list(title="Caudal (lps/km2)",range=c(-1,2),type = "log"),
             yaxis = list(title="NO3 (mg/s/km2)", range=c(-2,1),type = "log"),
             legend = list(x = 0.7, y = 0.01)
             #               showlegend = FALSE,
)


PT

fit = lm(log(NO3_mgs_km2)~log(flow_lps_km2), data=qlt.sim_23)
a = exp(fit$coefficients[1])
b = fit$coefficients[2]

q23fit = seq(min(qlt.sim_23$flow_lps_km2), max(qlt.sim_23$flow_lps_km2), length.out = 30)
NO3_23fit = a*q23fit^b


PT = plot_ly(data = qlt.sim_23, x=~flow_lps_km2,y=~NO3_mgs_km2, type='scatter', 
             mode='markers', name="Simulado",
             marker=list(color='gray', size=4))

PT = PT %>% add_markers(data=qlt_23,x = ~qobs, y = ~NO3_mgxs,
                        type = 'scatter',
                        mode = 'markers', 
                        marker = list(color = 'red', size=8),
                        name="Observado")

PT <- PT %>% add_lines(data=qlt_23,x = ~qobs, y = ~fit_NO3,
                       type = 'scatter',
                       mode = 'line', 
                       marker = list(color = 'red', size=0),
                       line = list(color = 'red', width = 2),
                       name="Regresión observados")

PT <- PT %>% add_lines(x= ~q23fit, y = ~NO3_23fit,
                       type = 'scatter',
                       mode = 'line', 
                       marker = list(color = 'black', size=1),
                       line = list(color = 'black', width = 2),
                       name="Regresión simulados")

PT <- layout(PT, xaxis = list(title="Caudal (lps/km2)",range=c(-1,2),type = "log"),
             yaxis = list(title="NO3 (mg/s/km2)", range=c(-2,1),type = "log"),
             legend = list(x = 0.6, y = 0.01)
             #               showlegend = FALSE,
)


PT

Histograma de parámetros

best_par <- readRDS("~/R_projects/SWAT_calibration_ETmodis/resultados/optim20200401_20/best_par.RDS")


par_name = names(par_lwr)[order(names(par_lwr))]

par(mfrow=c(12,2), mar=c(4,4,1,1))

for(i in 1:22){
p1 = par_name[i]
p2 = par_name[i+22]
hist(best_par[,p1], main=paste(p1),
     breaks = seq(par_lwr[p1], par_upr[p1], length.out = 30))

hist(best_par[,p2], main=paste(p2),
     breaks = seq(par_lwr[p2], par_upr[p2], length.out = 30))
}

Anexo

rm(list = ls())

# Librerias
library(SWATplusR)
library(stringr)
library(hydroGOF)
library(hydromad)
library(tibble)
library(hydroPSO)

# Caudales observados
obs.data <- readRDS("./data/obs.data.RDS")
obs.data = as.data.frame(obs.data)
obs.data$q23lag1 = c(obs.data$q23[2:nrow(obs.data)],NA)
obs.data$q36lag1 = c(obs.data$q36[2:nrow(obs.data)],NA)
obs.data$q41lag2 = c(obs.data$q41[3:nrow(obs.data)],NA,NA)

obs.data$PT_23_lag1 = c(obs.data$PT_23[2:nrow(obs.data)],NA)
obs.data$NO3_23_lag1 = c(obs.data$NO3_23[2:nrow(obs.data)],NA)
obs.data$PT_41_lag2 = c(lag(obs.data$PT_41)[3:nrow(obs.data)],NA,NA)
obs.data$NO3_41_lag2 = c(lag(obs.data$NO3_41)[3:nrow(obs.data)],NA,NA)


# SWAT project mensual
swat_project = "./SWAT_SL_20200401/"

# Funcion de optimización
swat_optim <- function(mypar, obs.data, project_dir, start_sim, end_sim,runrun_path, head_log) {
  
  path_run = paste0(runrun_path,"run_",Sys.getpid())
  
  par_set = tibble("r1_CN2::CN2.mgt|change = pctchg |
                   sub == c(1,8,10,14,15,16,17,18,19,20,22)" =
                     mypar[1]*100,
                   "r1_ESCO::ESCO.hru|change = absval |
                   sub == c(1,8,10,14,15,16,17,18,19,20,22)" =
                     mypar[2],
                   "r1_SOL_AWC::SOL_AWC.sol|change = pctchg |
                   sub == c(1,8,10,14,15,16,17,18,19,20,22)" =
                     mypar[3]*100,
                   "r1_ALPHA_BF::ALPHA_BF.gw|change = absval |
                   sub == c(1,8,10,14,15,16,17,18,19,20,22)" =
                     mypar[4],
                   "r1_GWQMN::GWQMN.gw|change = absval |
                   sub == c(1,8,10,14,15,16,17,18,19,20,22)" =
                     mypar[5],
                   "r1_GW_DELAY::GW_DELAY.gw|change = absval |
                   sub == c(1,8,10,14,15,16,17,18,19,20,22)" =
                     mypar[6],
                   "r1_REVAPMN::REVAPMN.gw|change = absval |
                   sub == c(1,8,10,14,15,16,17,18,19,20,22)" =
                     mypar[7],
                   "r1_GW_REVAP::GW_REVAP.gw|change = absval |
                   sub == c(1, 16, 18, 22,23)" =
                     mypar[8],
                   "r1_OV_N::OV_N.hru|change = absval |
                   sub == c(1,8,10,14,15,16,17,18,19,20,22)" =
                     mypar[9],
                   "r1_SLSUBBSN::SLSUBBSN.hru|change = pctchg |
                   sub == c(1,8,10,14,15,16,17,18,19,20,22)" =
                     mypar[10]*100,
                   
                   "r2_CN2::CN2.mgt|change = pctchg |
                   sub == c(2,3,4,5,6,7,9,11,12,13,21,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41)" =
                     mypar[11]*100,
                   "r2_ESCO::ESCO.hru|change = absval |
                   sub == c(2,3,4,5,6,7,9,11,12,13,21,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41)" =
                     mypar[12],
                   "r2_SOL_AWC::SOL_AWC.sol|change = pctchg |
                   sub == c(2,3,4,5,6,7,9,11,12,13,21,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41)" =
                     mypar[13]*100,
                   "r2_ALPHA_BF::ALPHA_BF.gw|change = absval |
                   sub == c(2,3,4,5,6,7,9,11,12,13,21,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41)" =
                     mypar[14],
                   "r2_GWQMN::GWQMN.gw|change = absval |
                   sub == c(2,3,4,5,6,7,9,11,12,13,21,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41)" =
                     mypar[15],
                   "r2_GW_DELAY::GW_DELAY.gw|change = absval |
                   sub == c(2,3,4,5,6,7,9,11,12,13,21,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41)" =
                     mypar[16],
                   "r2_REVAPMN::REVAPMN.gw|change = absval |
                   sub == c(2,3,4,5,6,7,9,11,12,13,21,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41)" =
                     mypar[17],
                   "r2_GW_REVAP::GW_REVAP.gw|change = absval |
                   sub == c(2,3,4,5,6,7,9,11,12,13,21,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41)" =
                     mypar[18],
                   "r2_OV_N::OV_N.hru|change = absval |
                   sub == c(2,3,4,5,6,7,9,11,12,13,21,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41)" =
                     mypar[19],
                   "r2_SLSUBBSN::SLSUBBSN.hru|change = pctchg |
                   sub == c(2,3,4,5,6,7,9,11,12,13,21,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41)" =
                     mypar[20]*100,
                   
                   "r1_HRU_SLP::HRU_SLP.hru|change= pctchg | 
                   sub == c(1,8,10,14,15,16,17,18,19,20,22)"= mypar[21]*100,
                   "r1_USLE_K::USLE_K.sol|change= pctchg | 
                   sub == c(1,8,10,14,15,16,17,18,19,20,22)" = mypar[22]*100,
                   "r1_USLE_P::USLE_P.mgt|change= pctchg | 
                   sub == c(1,8,10,14,15,16,17,18,19,20,22)"=mypar[23]*100,
                   "r1_P_UPDIS::P_UPDIS.bsn|change= absval | 
                   sub == c(1,8,10,14,15,16,17,18,19,20,22)"= mypar[24], 
                   "r1_FRT_SURFACE::FRT_SURFACE.mgt|change=absval | 
                   sub == c(1,8,10,14,15,16,17,18,19,20,22)" = mypar[25], 
                   "r1_FRT_KG::FRT_KG.mgt|change=pctchg | 
                   sub == c(1,8,10,14,15,16,17,18,19,20,22)" = mypar[26]*100,
                   "r1_NPERCO::NPERCO.bsn|change= absval | 
                   sub == c(1,8,10,14,15,16,17,18,19,20,22)"=mypar[27],
                   "r1_ERORGN::ERORGN.hru|change= absval | 
                   sub == c(1,8,10,14,15,16,17,18,19,20,22)"= mypar[28], 
                   
                   "r2_HRU_SLP::HRU_SLP.hru|change= pctchg | 
                   sub == c(2,3,4,5,6,7,9,11,12,13,21,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41)"= 
                     mypar[29]*100,
                   "r2_USLE_K::USLE_K.sol|change= pctchg | 
                   sub == c(2,3,4,5,6,7,9,11,12,13,21,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41)" = 
                     mypar[30]*100,
                   "r2_USLE_P::USLE_P.mgt|change= pctchg | 
                   sub == c(2,3,4,5,6,7,9,11,12,13,21,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41)"= 
                     mypar[31]*100,
                   "r2_P_UPDIS::P_UPDIS.bsn|change= absval | 
                   sub == c(2,3,4,5,6,7,9,11,12,13,21,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41)"= 
                     mypar[32], 
                   "r2_FRT_SURFACE::FRT_SURFACE.mgt|change=absval | 
                   sub == c(2,3,4,5,6,7,9,11,12,13,21,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41)"= 
                     mypar[33], 
                   "r2_FRT_KG::FRT_KG.mgt|change=pctchg | 
                   sub == c(2,3,4,5,6,7,9,11,12,13,21,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41)"= 
                     mypar[34]*100,
                   "r2_NPERCO::NPERCO.bsn|change= absval | 
                   sub == c(2,3,4,5,6,7,9,11,12,13,21,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41)"= 
                     mypar[35],
                   "r2_ERORGN::ERORGN.hru|change= absval | 
                   sub == c(2,3,4,5,6,7,9,11,12,13,21,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41)"= 
                     mypar[36], 
                   
                   "r1_PPERCO::PPERCO.bsn|change= absval | 
                   sub == c(1,8,10,14,15,16,17,18,19,20,22)"= mypar[37],                    
                   "r2_PPERCO::PPERCO.bsn|change= absval | 
                   sub == c(2,3,4,5,6,7,9,11,12,13,21,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41)"= 
                     mypar[38],
                   
                   "r1_CDN::CDN.bsn|change= absval | 
                   sub == c(1,8,10,14,15,16,17,18,19,20,22)"= mypar[39],                    
                   "r2_CDN::CDN.bsn|change= absval | 
                   sub == c(2,3,4,5,6,7,9,11,12,13,21,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41)"= 
                     mypar[40],
                   
                   "r1_CMN::CMN.bsn|change= absval | 
                   sub == c(1,8,10,14,15,16,17,18,19,20,22)"= mypar[41],                    
                   "r2_CMN::CMN.bsn|change= absval | 
                   sub == c(2,3,4,5,6,7,9,11,12,13,21,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41)"= 
                     mypar[42],
                   
                   "r1_PSP::PSP.bsn|change= absval | 
                   sub == c(1,8,10,14,15,16,17,18,19,20,22)"= mypar[43],                    
                   "r2_PSP::PSP.bsn|change= absval | 
                   sub == c(2,3,4,5,6,7,9,11,12,13,21,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41)"= 
                     mypar[44]
  ) 
  
  
  # run swat
  qsim = run_swat2012(project_path = project_dir,
                      output = list(FLOW = define_output(file = "rch",
                                                         variable = "FLOW_OUT",
                                                         unit = c(23,36,41)),
                                    PT = define_output(file = "rch",
                                                       variable = "TOT_Pkg",
                                                       unit = c(23,41)),
                                    NO3 = define_output(file = "rch",
                                                        variable = "NO3_OUT",
                                                        unit = c(23,41))),
                      parameter = par_set,
                      start_date = as.Date("2009-01-01"),
                      end_date = as.Date("2015-12-31"),
                      years_skip = 1,
                      run_path=path_run,
                      quiet = TRUE,
                      output_interval="d",
                      add_parameter = FALSE)
  
  # objective function
  nse23 = NSE(as.vector(qsim$FLOW_23), as.vector(obs.data$q23lag1))
  nse36 = NSE(as.vector(qsim$FLOW_36), as.vector(obs.data$q36lag1))
  nse41 = NSE(as.vector(qsim$FLOW_41), as.vector(obs.data$q41lag2))
  
  kge23 = KGE(as.vector(qsim$FLOW_23), as.vector(obs.data$q23lag1))
  kge36 = KGE(as.vector(qsim$FLOW_36), as.vector(obs.data$q36lag1))
  kge41 = KGE(as.vector(qsim$FLOW_41), as.vector(obs.data$q41lag2))
  
  pbias23 = pbias(as.vector(qsim$FLOW_23), as.vector(obs.data$q23lag1))
  pbias36 = pbias(as.vector(qsim$FLOW_36), as.vector(obs.data$q36lag1))
  pbias41 = pbias(as.vector(qsim$FLOW_41), as.vector(obs.data$q41lag2))
  
  
  # objective quality
  qlt_41 = data.frame(flow_lps_km2=1000*qsim$FLOW_41/4896,
                      PT_mgs_km2 = 1000000*qsim$PT_41/(24*60*60*4896),
                      NO3_mgs_km2 = 1000000*qsim$NO3_41/(24*60*60*4896),
                      qobs_lps_km2 = 1000*obs.data$q41lag2/4896,
                      PTobs_mgs_km2 = 0.18*(1000*obs.data$q41lag2/4896)^1.02,
                      NO3obs_mgs_km2 = 0.11*(1000*obs.data$q41lag2/4896)^1.41)
  
  qlt_41[qlt_41<=0] = NA
  qlt_41 = na.omit(qlt_41)
  
  qlt_23 = data.frame(flow_lps_km2=1000*qsim$FLOW_23/687,
                      PT_mgs_km2 = 1000000*qsim$PT_23/(24*60*60*687),
                      NO3_mgs_km2 = 1000000*qsim$NO3_23/(24*60*60*687),
                      qobs_lps_km2 = 1000*obs.data$q23lag1/687,
                      PTobs_mgs_km2 = 0.1*(1000*obs.data$q23lag1/687)^0.93,
                      NO3obs_mgs_km2 = 0.08*(1000*obs.data$q23lag1/687)^1.41)
  
  qlt_23[qlt_23<=0] = NA
  qlt_23 = na.omit(qlt_23)
  
  fit = lm(log(PT_mgs_km2)~log(flow_lps_km2), data=qlt_41)
  a = exp(fit$coefficients[1])
  b = fit$coefficients[2]
  qlt_41$PT_fit = a*(qlt_41$flow_lps_km2)^b
  
  fit = lm(log(NO3_mgs_km2)~log(flow_lps_km2), data=qlt_41)
  a = exp(fit$coefficients[1])
  b = fit$coefficients[2]
  qlt_41$NO3_fit = a*(qlt_41$flow_lps_km2)^b
  
  fit = lm(log(PT_mgs_km2)~log(flow_lps_km2), data=qlt_23)
  a = exp(fit$coefficients[1])
  b = fit$coefficients[2]
  qlt_23$PT_fit = a*(qlt_23$flow_lps_km2)^b
  
  fit = lm(log(NO3_mgs_km2)~log(flow_lps_km2), data=qlt_23)
  a = exp(fit$coefficients[1])
  b = fit$coefficients[2]
  qlt_23$NO3_fit = a*(qlt_23$flow_lps_km2)^b
  
  kge_PT_41 = KGE(qlt_41$PT_fit, qlt_41$PTobs_mgs_km2)
  kge_NO3_41 = KGE(qlt_41$NO3_fit, qlt_41$NO3obs_mgs_km2)
  kge_PT_23 = KGE(qlt_23$PT_fit, qlt_23$PTobs_mgs_km2)
  kge_NO3_23 = KGE(qlt_23$NO3_fit, qlt_23$NO3obs_mgs_km2)
  
  pbias_PT_23 = pbias(qlt_23$PT_fit, qlt_23$PTobs_mgs_km2)
  pbias_PT_41 = pbias(qlt_41$PT_fit, qlt_41$PTobs_mgs_km2)
  
  pbias_NO3_23 = pbias(qlt_23$NO3_fit, qlt_23$NO3obs_mgs_km2)
  pbias_NO3_41 = pbias(qlt_41$NO3_fit, qlt_41$NO3obs_mgs_km2)
  
  
  FO23 = 0.5*kge23 + 0.1*kge_NO3_23 + 0.1*kge_PT_23 + 0.15*(1-abs(pbias_PT_23)/100) + 
    0.15*(1-abs(pbias_NO3_23)/100)
  FO41 = 0.5*kge41 + 0.1*kge_NO3_41 + 0.1*kge_PT_41 + 0.15*(1-abs(pbias_PT_41)/100) + 
    0.15*(1-abs(pbias_NO3_41)/100)
  FO = 0.6*FO41 + 0.4*FO23 
  
  # logfile
  vout = data.frame(FO,kge23,kge36,kge41,
                    nse23,nse36,nse41,
                    pbias23,pbias36,pbias41,
                    kge_PT_23, kge_PT_41, kge_NO3_41, kge_NO3_23,
                    pbias_PT_23, pbias_PT_41, pbias_NO3_23, pbias_NO3_41,
                    t(mypar))
  colnames(vout) = head_log
  log.file = paste0(runrun_path,"log_",Sys.getpid(),".txt")
  
  if(!any(list.files(runrun_path)==paste0("log_",Sys.getpid(),".txt"))){
    write.table(t(head_log), log.file, col.names = FALSE, row.names = FALSE, quote = FALSE)
  }
  
  ff = read.table(log.file, col.names = head_log,colClasses = "character")
  ff = rbind(ff,vout)
  write.table(ff, log.file, col.names = FALSE, row.names = FALSE, quote = FALSE)
  
  #output
  return(FO)
}

r1_CN2 = c(-0.4, 0)   ###
r2_CN2 = c(-0.2,0.2)
r1_ESCO = c(0.1,0.4)
r2_ESCO = c(0.35,0.45)
r1_SOL_AWC = c(0,0.2)
r2_SOL_AWC = c(-0.15,0)
r1_ALPHA_BF = c(0.4,0.8)  
r2_ALPHA_BF = c(0.4,0.8)  
r1_GWQMN = c(800,6000)
r2_GWQMN = c(800,1500)
r1_GW_DELAY = c(20,35)
r2_GW_DELAY = c(25,45)
r1_REVAPMN = c(325,375)
r2_REVAPMN = c(610,710)
r1_GW_REVAP = c(0.02, 0.1)   ###
r2_GW_REVAP = c(0.03, 0.055)
r1_OV_N = c(0.4, 0.6)
r2_OV_N = c(0.4, 0.6)
r1_SLSUBBSN = c(-0.2, 0.05) ###
r2_SLSUBBSN = c(0.2, 0.4)

r1_HRU_SLP = c(-0.1,0.1)
r2_HRU_SLP = c(-0.1,0.1)
r1_USLE_K = c(-0.15,0.15)
r2_USLE_K = c(0,0.2)
r1_USLE_P = c(-0.15,0.15)
r2_USLE_P = c(-0.15,0.15)
r1_P_UPDIS = c(20, 60)
r2_P_UPDIS = c(40, 90)
r1_FRT_SURFACE = c(0.1,0.6)
r2_FRT_SURFACE = c(0.1,0.6)
r1_FRT_KG = c(-0.45,-0.15)
r2_FRT_KG = c(-0.4,0.1)
r1_NPERCO = c(0.5,0.8)
r2_NPERCO = c(0.7,1)
r1_ERORGN = c(3,4)
r2_ERORGN = c(2,4)

r1_PPERCO = c(10,17.5)
r2_PPERCO = c(10,17.5) #bsn absval
r1_CDN = c(0,3) #bsn absval
r2_CDN = c(0,3) #bsn absval
r1_CMN = c(0,0.03) #bsn absval
r2_CMN = c(0,0.03) #bsn absval
r1_PSP = c(0.01,0.8) #bsn absval
r2_PSP = c(0.01,0.8)

#           CN2   ESCO  SOL_AWC ALPHA_BF  GWQMN GW_DELAY  REVAPMN GW_REVAP  OV_N  SLSUBBSN
par_lwr = c(r1_CN2[1], r1_ESCO[1], r1_SOL_AWC[1],r1_ALPHA_BF[1], r1_GWQMN[1], r1_GW_DELAY[1], 
            r1_REVAPMN[1], r1_GW_REVAP[1], r1_OV_N[1], r1_SLSUBBSN[1],
            r2_CN2[1], r2_ESCO[1], r2_SOL_AWC[1],r2_ALPHA_BF[1], r2_GWQMN[1], r2_GW_DELAY[1], 
            r2_REVAPMN[1], r2_GW_REVAP[1], r2_OV_N[1], r2_SLSUBBSN[1],
            r1_HRU_SLP[1], r1_USLE_K[1], r1_USLE_P[1], r1_P_UPDIS[1], r1_FRT_SURFACE[1], 
            r1_FRT_KG[1], r1_NPERCO[1], r1_ERORGN[1],
            r2_HRU_SLP[1], r2_USLE_K[1], r2_USLE_P[1], r2_P_UPDIS[1], r2_FRT_SURFACE[1], 
            r2_FRT_KG[1], r2_NPERCO[1], r2_ERORGN[1],
            r1_PPERCO[1], r1_CDN[1], r1_CMN[1], r1_PSP[1],
            r2_PPERCO[1], r2_CDN[1], r2_CMN[1], r2_PSP[1])

par_upr = c(r1_CN2[2], r1_ESCO[2], r1_SOL_AWC[2],r1_ALPHA_BF[2], r1_GWQMN[2], r1_GW_DELAY[2], 
            r1_REVAPMN[2], r1_GW_REVAP[2], r1_OV_N[2], r1_SLSUBBSN[2],
            r2_CN2[2], r2_ESCO[2], r2_SOL_AWC[2],r1_ALPHA_BF[2], r2_GWQMN[2], r2_GW_DELAY[2], 
            r2_REVAPMN[2], r2_GW_REVAP[2], r2_OV_N[2], r2_SLSUBBSN[2],
            r1_HRU_SLP[2], r1_USLE_K[2], r1_USLE_P[2], r1_P_UPDIS[2], r1_FRT_SURFACE[2],
            r1_FRT_KG[2], r1_NPERCO[2], r1_ERORGN[2],
            r2_HRU_SLP[2], r2_USLE_K[2], r2_USLE_P[2], r2_P_UPDIS[2], r2_FRT_SURFACE[2], 
            r2_FRT_KG[2], r2_NPERCO[2], r2_ERORGN[2],
            r1_PPERCO[2], r1_CDN[2], r1_CMN[2], r1_PSP[2],
            r2_PPERCO[2], r2_CDN[2], r2_CMN[2], r2_PSP[2])

par_init = c(-0.198073840010168,0.286267953999935,0.00746275760208839,
             0.499157443870258,3500.6256700673,33.048081707854,339.59120477244,
             0.032221701353556,0.514936793955879, -0.00241425477305834,0.0219494794735184,
             0.41861314804868,-0.113031946972844, 0.69116116427558,1456.88879168975,39.7176103641318,
             660.900428427947,0.0350548241647086, 0.559651736719884,0.380954113796417,
             0.0156559405006919,0.0968484356612427,-0.0331742106779189,40.6194318342368,
             0.386746499579996,-0.34870424019712,0.646162500043787,3.65240741158593,0.0113505467880799,
             0.139763966103355,0.0545959533479474,50.2166504226717,0.337105246566447,-0.23731772825072,
             0.886945397111459,3.03117861658487,13, 1.1,0.01,0.2, 13, 1.1,0.01,0.2)

# Optimización
head_log = c("FO","KGE23","KGE36", "KGE41",
             "NSE23","NSE36", "NSE41",
             "pbias23","pbias36", "pbias41",
             "kge_PT_23", "kge_PT_41", "kge_NO3_41", "kge_NO3_23",
             "pbias_PT_23", "pbias_PT_41", "pbias_NO3_23", "pbias_NO3_41",
             "r1_CN2", "r1_ESCO", "r1_SOL_AWC","r1_ALPHA_BF","r1_GWQMN", 
             "r1_GW_DELAY","r1_REVAPMN","r1_GW_REVAP","r1_OV_N","r1_SLSUBBSN",
             "r2_CN2", "r2_ESCO", "r2_SOL_AWC","r2_ALPHA_BF","r2_GWQMN", 
             "r2_GW_DELAY","r2_REVAPMN","r2_GW_REVAP","r2_OV_N","r2_SLSUBBSN",
             "r1_HRU_SLP","r1_USLE_K","r1_USLE_P","r1_P_UPDIS","r1_FRT_SURFACE",
             "r1_FRT_KG","r1_NPERCO","r1_ERORGN",
             "r2_HRU_SLP","r2_USLE_K","r2_USLE_P","r2_P_UPDIS","r2_FRT_SURFACE",
             "r2_FRT_KG","r2_NPERCO","r2_ERORGN",
             "r1_PPERCO", "r1_CDN", "r1_CMN", "r1_PSP",
             "r2_PPERCO", "r2_CDN", "r2_CMN", "r2_PSP")


fr = list.files("./run_optim_mensual/", pattern = "*.txt", full.names=TRUE)
file.remove(fr)
opt_pso = hydroPSO(par=par_init, fn = swat_optim,lower=par_lwr, upper=par_upr, 
                   control = list(MinMax = "max",
                                  maxit=500,
                                  reltol=0.00001,
                                  parallel="parallelWin", 
                                  par.nnodes=12,
                                  par.pkgs = list("hydromad", 
                                                  "SWATplusR",
                                                  "stringr",
                                                  "hydroGOF",
                                                  "tibble"),
                                  write2disk=FALSE,
                                  REPORT=1,
                                  npart=48,
                                  normalise=TRUE),
                   obs.data = obs.data, 
                   project_dir = swat_project, 
                   runrun_path = "./run_optim_mensual/",
                   head_log = head_log
)